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ABSTRACT 


The  present  paper  studies  -the-  linear  complementarity  problem  v 

of  finding  vectors  x and  y in  such  that  c + Dx  + y > 0,  J 

T T 

b - x > 0 and  x (c  + Dx  + y)  = y (b  - x)  - 0 where  D is  a 

^ vv  h t cK\x 

Z-matrix  and  b > 0.  Complementarity  problems  of  this  nature  aris^j 

for  example^  fro®  "the  minimization  of  certain  quadratic  functions  subject 
to  upper  and  lower  bounds  on  the  variables.  Two  least-element  character- 
izations of  solutions  to  the  above  linear  complementarity  problem  are 
established  first.  Next,  a new  and  direct  method  to  solve  this  class 
of  problems,  which  depends  on  the  idea  of  /least- element  solution^is 


presented.  Finally,  applications  and  computational  experience  with 


i 


1.  INTRODUCTION 

The  present  paper  is  concerned  with  the  linear  complementarity 
problem  of  finding  vectors  x,  y £ R^  such  that 

u =c  + Dx  + y>0 

(l.l)  v = b - x > 0 

T T. 
u x = v y = 0 

where  b £ R^>  c £ Rn  and  D £ RnXn.  We  denote  problem  (l.l)  by  the 
triple  (b,  c,D).  If  D is  symmetric,  which  incidentally,  is  not  assumed 


in  this  paper,  then  (l.l)  is  precisely  the  Kuhn-Tucker  optimality 
conditions  for  the  quadratic  program  of  finding  a vector  x £ R°  to 

T IT 

(1.2)  minimize  c x + i x Dx  subject  to  x < b . 

It  is  clear  that  any  quadratic  program  of  minimizing  a quadratic  function 
subject  to  upper  and  lower  bounds  on  the  variables  can  be  cast  in  the 
form  (1.2). 

The  theory  and  applications  of  the  linear  complementarity 
problems  with  Z-matrices  (to  be  defined  in  the  next  section),  which, 
with  the  assumption  of  symmetry  on  the  matrices,  correspond  to  the 
minimization  of  certain  quadratic  functions  subject  to  lower  bounds  on 
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the  variables,  have  received  much  attention  in  the  literature  [2], 
[7]>  [8],  [9],  [10],  [18],  [19],  [22].  Recently,  many  applications 
of  both  problems  (1.1)  and  (1.2)  have  appeared  in  various  contexts; 


to  mention  a few,  the  unilateral  Dirichlet  problem  with  two  constraints 
[21]  leads  to  problem  (l.l)  where  D is  a Z-matrix;  the  taut  string 
problem  [23]  which  has  its  own  applications  in  inventory  theory  and 
statistics,  is  a special  case  of  (1.2)  where  the  matrix  D is 


Minkowski  and  tridiagonal;  and  Cheng's  salary  administration  model  [3] 

gives  rise  to  a problem  of  the  form  (1.2)  where  the  matrix  U = nl  - e eT 

n n n 

with  In  the  identity  matrix  of  order  n and  e^  = (l, ...,1)  £ Rn. 

In  all  these  instances  (and  many  others)  the  matrix  D 3 Z. 


Several  recent  papers  ([7],  [8],  [9],  [12],  [16],  [22])  have 
demonstrated  that  many  linear  complementarity  problems  have  solutions 


which  are  ’’least  element^  of  subsets  of  Euclidean  space.  In  an  earlier 
paper  [7],  R.  W.  Cottle  and  the  author  summarized  this  least-element 
aspect  of  the  solutions  for  various  classes  of  linear  complementarity 
problems.  We  focused  on  the  class  0 of  square  matrices  which  was 
introduced  by  Mangasarian  [14]  in  formulating  the  linear  complementarity 
problems  as  linear  programs  and  which  consists  of  square  matrices 
satisfying  the  two  conditions:  (i)  MX  = Y and  (ii)  rTX  + sTY  > 0 

for  some  r,  s € R^  and  where  X and  Y are  Z-matrices.  We 
demonstrated  that  linear  complementarity  problems  with  matrices  in  o 
have  solutions  which  can  be  obtained  as  least  elements  of  polyhedral 
sets  and  that  the  class  C contains  all  the  other  classes  of  matrices 
considered  in  [14],  [15]  and  in  particular,  all  those  previously  known 
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classes  of  matrices  ([8],  [9],  [12],  [22])  for  which  this  solution 


characterization  holds.  Despite  the  fact  that  0 contains  many 


interesting  matrices,  e.g.  those  enumerated  in  the  last  table  in  [15] 


not  too  many  realizations  of  the  linear  complementarity  problem  which 


yield  matrices  belonging  to  C have  been  seen.  Of  course,  as  mentioned 
earlier,  complementarity  problems  with  Z-matrices  (which  clearly  belong 
to  C)  have  many  applications,  e.g.  see  [18].  Later  in  the  paper,  we 
shall  show  that  the  matrix  ( ^ *),  which  appears  in  the  problem 

(b, c,D),  belongs  to  C if  and  only  if  D is  a Minkowski  matrix. 

In  many  potential  applications  of  problem  (l.l)  to  obtain 


numerical  solutions  of  partial  differential  equations  in  their  dis 


cretized  form,  the  matrix  D is  very  large,  sparse  and  structured,  in 


addition  to  being  of  class  Z.  See  [6].  Special  (iterative)  methods 


under  various  additional  assumptions  on  the  matrix  D have  been  proposed 
and  implemented  [1],  [6],  [15],  [17],  [21].  Direct  methods  like  the 
principal  pivoting  scheme  [4]  and  Lemke's  almost  complementarity  pivot- 


one  thing,  they  do  not  take  advantage  of  the  special  structure  that 


D possesses,  thus  creating  storage  difficulties  when  handling  large- 
scale  problems.  Cheng  [5]  and  Veinott  [23]  have  described  special  (direct) 


methods  for  solving  the  salary  administration  problem  and  the  taut 


string  problem  respectively.  Unfortunately,  their  methods  are  not 


applicable  to  the  general  problem. 

Our  purposes  in  this  paper  are:  (i)  to  establish  two  least- 

element  aspects  of  the  linear  complementarity  problem  (b,c,D); 


(ii)  to  present  a new  (direct)  method  for  solving  large-scale  linear 
complementarity  problems  of  this  class;  and  (iii)  to  report  our 
computational  experience  in  solving  some  realizations  of  this  class 
of  problems  using  this  new  method.  The  plan  of  the  paper  is  as  follows. 
In  the  next  section,  we  review  some  basic  terminology,  fix  our  notations 
and  present  a result  in  lattice  theory  that  will  be  used  in  later 
development.  In  Section  5,  we  develop  the  least-element  aspects  of 
the  problem.  In  Section  k,  we  present  our  proposed  method  and  discuss 
some  of  its  refinements  when  it  is  applied  to  problems  having  further 
structure.  In  the  fifth  and  final  section  of  this  paper,  we  report 
our  computational  experience  with  the  method. 
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2.  NOTATIONS,  DEFINITIONS  AMD  BASIC  RESULTS 


We  denote  by  r"  the  nonnegative  orthant  of  Euclidean  n-space. 

By  Rnxm,  we  denote  the  space  of  real  nxm  matrices.  We  use  I to 
denote  the  identity  matrix  of  order  n and  eR  to  denote  the  summation 
vector  (l,  € Rn.  The  real  matrix  A £ RnXn  is  said  to  be  a 

Z -matrix  (P-matrix)  if  it  has  nonpositive  off-diagonal  entries  (positive 
principal  minors ) . We  shall  call  a matrix  A £ Rn  x n a K-matrix  (or  a 
Minkowski  matrix)  if  it  is  both  a Z-  and  a P-matrix  simultaneously. 

The  classes  of  all  real  Z-,  P-  and  K-matrices  will  be  denoted  by  Z, 

P and  K respectively.  It  is  obvious  that  principal  submatrices  of 
Z-,  P-  and  K-matrices  are  themselves  Z-,  P-  and  K-matrices  respectively. 
Proofs  of  the  following  characterizations  of  P-  and  K-matrices  can  be 
found  in  Fiedler  and  Ptak  [ 11] . 
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Proposition  2.1.  (i)  A matrix  M is  a P-matrix  if  and  only  if  to 

every  vector  x ^ 0,  there  exists  an  index  k such  that  x^Mx^X). 
(ii)  Let  M £ Z.  Then  M £ K if  and  only  if  there  exists  a vector 
x > 0 such  that  Mx  > 0. 

For  a vector  q £ Rn  and  a matrix  M £ Rn  X n,  the  linear  com- 
plementarity problem,  (q,M)  is  that  of  finding  x £ R^  such  that 

q + Mx  > 0 and  x^(q  + Mx)  = 0 . 

Note  that  the  problem  (b,c,D)  is  precisely  the  linear  complementarity 


problem  (q,M)  with  q = (£)  and  M = (Jj.  ^).  The  problem  (q,M) 
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is  said  to  be  feasible  if  there  exists  a vector  x € R 


such  that 


vectors  is  called  the  feasible  set.  It  has  been  shown  (see  Samelson 


et  al.  [20])  that  (q,M)  has  a unique  solution  for  every  q £ Rn  if 


is  a principal  submatrix  of  M.  Similarly,  for  a vector  q € R 


For  the  sake  of  completeness,  we  briefly  review  a few  concepts 


in  lattice  theory  and  state  a theorem  which  is  essential  in  the  least 


element  study  of  the  problem  under  consideration.  The  following  dis- 
cussion and  a proof  of  the  theorem  can  be  found  in  Veinott  [2^]  and 
Pang  [16], 

A subset  S of  Rn  is  called  a meet  semi-sublattice  (of  Rn) 
if  for  every  x,  y € S,  the  vector  z = min(x,y)  defined  by 
z = min(x.,y, ) for  each  i,  also  belongs  to  S.  The  subset  S is 


bounded  below  if  there  exists  a vector  x'  € R such  that  x > x 


for  every  x € S.  An  element  x £ S is  a least  element  of  S if 
x > x for  every  x £ S.  It  is  clear  that  a least  element,  if  it 
exists,  must  be  unique. 

Theorem  2.2.  Let  S be  a nonempty  and  closed  meet  semi -sublattice 
of  Rn.  Suppose  that  S is  also  bounded  below.  Then  S has  a least 
element. 


Throughout  this  paper,  we  assume  D £ Z and  b > 0.  It  has 
been  proved  that  the  linear  complementarity  problem  with  a Z-matrix 
has  a solution  which  is  the  least  element  of  the  feasible  set,  provided 
that  the  latter  is  nonempty.  In  fact,  this  assertion  follows  immediately 
from  Theorem  2.2  by  noting  that  the  feasible  set  is  a meet  semi -sub- 
lattice of  Rn  satisfying  the  conditions  in  the  theorem.  It  is  then 
not  difficult  to  verify  that  the  least  element  solves  the  linear  com- 
plementarity problem.  We  note  that  even  though  D £ Z,  the  matrix 
( j q)  which  appears  in  the  problem  (b,c,D)  does  not  belong  to  Z. 
Therefore  the  above  assertion  does  not  apply.  Furthermore,  the  feasible 
set  of  the  problem,  which  is  defined  as  the  set 

{ (x)  € R^  x r“  : x < b and  c + Dx  + y > 0} 

y 

2n 

might  not  itself  be  a meet  semi -sublattice  of  R , as  easily  con- 
structed examples  will  show.  Therefore,  Theorem  2.2  cannot  be  applied 
to  the  above  set.  In  this  section,  we  develop  two  least  element  aspects 
of  the  problem.  Both  of  them  will  show  that  the  problem  possesses  a 
solution  which  can  be  obtained  from  the  least  element  of  subsets  of 
Euclidean  spaces,  one  of  which  is  in  R^  (the  x-space)  while  the  other 
in  r"  X R^  (the  space  of  the  feasible  set). 

We  start  our  least  element  study  of  the  problem  by  observing 
that,  even  without  the  assumption  D £ Z,  a vector  x £ R^  satisfying 
x < b clearly  determines  a (unique)  vector  y £ Rn  such  that  (x,y) 
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solves  the  problem  provided  the  conditions  below  are  satisfied  for 
each  i = 1, . . . , n: 

(i)  x^  = 0 =*(c  + Dx)^  > 0 

(ii)  b > xi>0  =*  (c  + Dx)i  = 0 

(iii)  ^ =>(c  + Dx)i  < 0 . 

With  an  abuse  of  language,  we  also  say  that  a vector  x £ R^  satisfying 
x < b and  conditions  (i)-(iii)  is  a solution  with  the  understanding 
that  the  y vector  exists  such  that  (x,y)  is  a real  solution. 

We  now  describe  a meet  semi-sublattice  of  Rn  having  a least 
element  which  solves  the  problem.  In  the  next  section,  we  will  present 
an  algorithm  which  actually  computes  this  least  element. 

Theorem  3.1.  Let  D £ Z and 
n 

S = H {x  £ R^:x  < b;  x^  < b^  =*  (c  + Dx)^  > 0} 
i=l 

Then  S is  a meet  semi-sublattice  of  Rn  and  has  a least  element. 
Moreover,  this  least  element  is  a solution  for  the  problem  (b,c,D). 

Proof;  Since  intersections  of  meet  semi-sublattices  are  meet  semi- 
sublattices, it  suffices  to  show  that  each  set 

S^  = {x  £ R^:x  < b;  x^  < b^  (c  + Dx) ^ > 0) 

is  one  such  for  i = 1, ,n.  Let  x,  x'  £ S^  and  x"  = min(x,x’). 

Clearly,  0 < x"  < b.  Suppose  x^  < ^ and  say,  x'^  = x^.  Then 
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<c  * Dx">i  ■ Ci  + diiXi  + L dijXJ 


oA 

E 

dA 


> o,  + d, ,x,  + E d. .x.  > 0 . 
- i ii  i ij  j - 


Thus  x"  £ S^.  Therefore,  S is  indeed  a meet  semi-sublattice  of  Rn. 
It  is  clearly  closed  and  bounded  below  by  0:  it  is  nonempty  because 
b £ S.  Hence  by  Theorem  2.2,  it  has  a least  element,  say  x.  It 
remains  to  verify  that  conditions  (ii)  and  (iii)  are  satisfied  for  x. 

We  omit  these  proofs  because  they  are  similar  to  the  one  given  in 
Lemma  3.10  in  Cottle  and  Pang  [71*  □ 


Theorem  3.1  above  shows  that  when  D £ Z,  the  problem  (b,c,D) 
has  a solution.  The  next  proposition  is  concerned  with  the  uniqueness 
of  the  solution. 

Proposition  3.2.  Suppose  D £ K.  Then  the  problem  (b,c,D)  has  a 
unique  solution  (x,y)  where  x is  the  least  element  of  S . 

Proof:  It  suffices  to  establish  the  uniqueness  part.  Let  (x,y)  be 
another  solution.  Then  for  i = 1, ...,n,  we  have 

(x  - x^Cu  - u^  = (x  - xj^DCx  - x)  + (y  - y)Ji 

= (x  - x)i  (D(x  - x))i  + (x  - x)i(y  - y)i  . 

Furthermore, 
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= viyi  + viyi  > 0 . 

Thus  (x  - x)i(D(x  - x))i  < 0 for  each  i.  Hence  it  follows  from 

Proposition  2.1  that  x = x.  Clearly,  y = y.  This  completes  the  proof. 

□ 

Corollary  3.3.  Suppose  D t K.  Let  x and  x he  the  (unique) 
solutions  of  the  problems  (b,c,D)  and  (c,D)  respectively.  Then  x<x. 

Proof:  Let  x = min(x,b).  Then  x S St  for  clearly  0 < x < b. 

If  x.^  < b^  then  x,^  = x^,  so  that 
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Therefore,  the  meet  semi-sublattice  S and  the  feasible  set  belong  to 


two  different  spaces.  Here,  our  objective  is  to  establish  that  the 


feasible  set  itself  contains  "some"  least  element,  i.e.  a least  element 


under  a different  partial  ordering,  that  solves  the  problem.  To  achieve 


this,  we  prove  the  theorem  below 


Let  D £ Z.  Then  D £ K 


Proof » Necessity.  If 


then  clearly  MX  = Y.  It  is  obvious  that  X £ K and  Y £ K.  Thus 


Then  there  exist 


Sufficiency 


r,  s such  that  MX  = Y 


and  rX+sY>0.  We  may  write 


are  nonnegative.  Then,  we  have 


Since  Y is  a Z-matrix,  it  follows  that  and  are  nonnegative 

T T T T T T 

diagonal  matrices.  We  may  write  r = (ri>r2)  and  s = (s*,  s*).  By 
an  easy  calculation,  we  obtain 


rTX  + sTY  = (rj,r£) 


X11  *12  \ „ m /“  11-^1  -«12  + ^2 


-X  X 

21  22 


/ T Tv 

+ (s^Sg) 


- (rixn  * #21  * *Xi  ' six2i  - 4xii> 


- rlX12  + f2X22  - =>12  + S1X22  + S2X12> 


Thus,  we  have 


(3.D 


rIXll  ' r2X21  + »Xl  - S1X21  - S2X11  > 0 


(3.2)  “rlX12  + #22  - #*12  + S1X22  + S2X02  > 0 


Expression  (3.1)  implies 


<rl  * #>  - 4>XU  > #21  + S1X21  ^ 0 ' 


Thus, 


T T T 
rL  + sjD  - sj  > 0 


Expression  (3.2)  implies 


fp  fp  fp  fp  fp 

(r"  + s})X22  > + s^D  - s^X^  > 0 


Since  *22  £ Z,  Proposition  2.l(ii)  implies  X22  € K.  The  fact  that 
Y £ Z implies  DX12  > X22.  Since  D £ Z and  X^  is  a nonnegative 
diagonal  matrix,  it  follows  that  DXly  £ Z.  Proposition  2.l(ii)  then 
implies  DX  ~ £ K and  thus  D £ K.  This  completes  the  proof.  t 


The  hypothesis  D £ Z is  necessary  in  order  for  the  sufficiency 


part  of  the  Theorem  3.4  to  be  true.  In  the  following  we  give  an  example 
of  a matrix  D such  that  D ^ Z but  ( ^ ?)  £ C . 


).  Obviously  D ^ Z.  Let  M = ( 


Choose 


Then  it  is  trivial  to  verify  MX  = Y.  Clearly  X,  Y £ Z.  If 
rT  = (4, 2, 3,1),  and  sT  = 0,  then  rTX  + sTY  > 0.  Therefore  M £ C 
Using  a result  (Theorem  3.H)  established  in  [7],  we  deduce 
that  when  D £ K,  the  (unique)  solution  of  the  problem  (b, c,D)  is 


given  by  the  vector 


is  the  least  element  of  the  polyhedral  set 


where 


Furthermore,  it  was  shown  in  the  same  reference  that  the  solution 


can  he  obtained  by  solving  the  linear  program  of  finding  vectors 


which 


minimize 


subject  to  c + Dx  + y > 0 


where  p and  q are  vectors  satisfying 


obtained  by  first  setting  p *=  r and  then  solving  D a - p + s for  q. 

The  discussion  above  shows  that  when  D £ K,  it  is  possible  to 
solve  the  linear  complementarity  problem  (b, c,D)  via  the  linear  program 
(3.3).  The  latter  can  be  solved  by  the  Simplex  Method  (with  the  upper 
bounding  technique)  or  the  iterative  (relaxation)  methods  for  systems 


of  linear  inequalities.  On  the  one  hand,  the  Simplex  Method  being 


very  general,  encounters  storage  difficulties  with  large  scale  problems 
of  this  nature.  On  the  other  hand,  our  computational  experience  with 


one  particular  iterative  method  is  quite  discouraging,  even  for  problems 
of  very  small  size  (see  [7]).  Although  further  analysis  and  clever 


modification  might  bring  about  improvements  in  the  computational 


performance  of  these  methods,  but  we  believe  that  another  approach  is 
preferable  for  solving  the  problem  (b,c,Dl.  In  the  next  section,  we 


propose  a new  and  efficient  algorithm  for  solving  large-scale  linear 


complementarity  problems  of  this  class.  It  is  a direct  method  which 


exploits  the  structure  of  the  problem, 


k.  A FAST  ALGORITHM  FOB  LARGE-SCALE  PROBLEMS 


Methods  for  solving  the  linear  complementarity  problem  (b, c,D) 
under  various  assumptions  (e.g.  symmetry  and  positive  definiteness) 
on  the  matrix  D have  been  proposed  and  implemented.  See  [6]  for  a 
survey.  These  include  both  iterative  procedures  [l] , [13],  [17],  [21] 
for  the  general  problem  and  direct  methods  for  some  of  its  particular 
cases  [3],  [23].  Typically,  in  applications  of  this  class  of  problems 
to  partial  differential  equations,  the  matrices  D are  very  large  and 
sparse.  The  kinds  of  direct  (pivoting)  methods  that  are  normally  used 
to  solve  linear  complementarity  problems  are  undesirable  because  they 
do  not  take  advantage  of  the  special  and  sparse  structure  of  the 
matrices,  thus  causing  storage  difficulties  when  handling  large-scale 
problems.  These  features  (special  structure  and  sparsity)  are  very 
important  factors  motivating  the  design  of  a special  efficient  algorithm. 
Our  purpose,  in  this  section  is  to  present  a new  (direct)  algorithm 
for  large-scale  linear  complementarity  problems  of  this  nature.  In 
the  first  part  of  the  section,  we  formulate  the  algorithm  in  its  general 
form.  In  the  latter  part  of  the  section,  we  refine  the  algorithm  to 
solve  the  subclass  of  problems  having  tridiagonal  matrices  D and  a 
class  of  variance-minimization  problems. 

Because  of  the  similarity  of  our  algorithm  and  Chandrasekaran' s 
algorithm  [2]  for  linear  complementarity  problems  with  Z-matrices,  we 
begin  by  reviewing  the  statement  of  the  latter  algorithm.  It  can  be 
shown  that  the  solution  obtained  by  this  algorithm  is  the  least  element 
of  the  feasible  set  [23]. 


17 


Algorithm  I:  Chandrasekaran' s Algorithm  for  (q,M)  with  M Z. 


Step  0.  Let  k = 0 and  x ^ = 0. 

Step  1.  Let  = [i;(q  + Mx^k))i  < 0}  and  J U ' = {1, . . . ,n  }\l (k  If 

1^  = (jS  stop.  A solution  is  x^.  Otherwise  continue. 

Step  2.  Let  I = I^k)  and  J = J^k).  Solve  Mnx^k+1)  = -q  If 

this  system  of  equations  does  not  possess  a solution,  stop. 

The  problem  (q,M)  is  infeasible.  Otherwise  set  xlk+1^  = 0. 

J 

Replace  k by  k+1  and  go  to  Step  1. 

In  what  follows,  we  propose  a slight  modification  of  the  above 
algorithm  for  the  case  M € K. 

Algorithm  II:  Modified  Chandrasekaran  Algorithm  for  (q,M)  with  M € K. 


Step  0. 

Let  k = 0 and  x^  = 0. 

Step  1. 

If  (q  + Mx^)  > 0,  stop.  The  solution  is  x^. 

Otherwise, 

let  I^k)  = {i:(q  + Mx^k))i  < 0)  and  J^={1,.. 

Continue. 

Step  2. 

Let  I = 1^,  and  J = J^.  Solve  M^x^4^  = - 
x^k+!)  _ Replace  k by  k+1  and  go  to  Step  1. 

q^  and  set 

The  two  algorithms  differ  in  the  definition  of  the  set  I^k' . 

In  the  former,  the  set  I^k-  is  defined  as  { i : ( q + Mx^)^  < 0} 
which,  in  general  is  more  restrictive  than  the  one  defined  in  the  latter. 


... 
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In  order  to  explain  the  motivation  for  our  modification,  we  recall  that 
an  essential  purpose  of  Chandrasekaran' s algorithm  is  to  identify  the 
set  I^  of  all  those  indices  i for  which  (q  + Mx)  ^ = 0 where  x 
is  a solution  (provided  that  it  exists)  to  the  problem.  This  purpose 
is  achieved  by  using  Chandrasekaran' s observation  (q  + Mx'  ')^<0=»£^>0 

which,  together  with  the  complementarity  condition,  implies  (q  + Mx)i  = 0. 

(k) 

Thus  the  sequence  of  sets  I forms  a successive  approximation  to 

the  desired  set  I^.  Recognizing  this  fact,  we  see  that  the  algorithm 

(k) 

may  be  improved  if  the  sets  Iv  can  be  made  to  "converge"  faster 

to  1^,  or  in  other  words,  if  fewer  systems  of  equations  have  to  be 

solved  in  Step  2.  We  therefore  propose  to  modify  the  algorithm  slightly 

by  observing  that  this  stronger  implication  holds,  namely, 

(q  + Mx^)^  < 0 =>  (q  + Mx)^  = 0,  and  redefine  the  set  1^  as  in 

(k)  - 

the  second  algorithm.  Note  that  the  implication  (q  + Mxv  <0=>x^>0 

(ir\ 

does  not  necessarily  hold.  The  new  definition  of  I ' is  intended  to 
include  as  many  indices  in  as  possible.  It  has  the  potential 

advantage  of  aggregating  several  systems  of  equations  into  a single 
one,  therefore  reducing  the  number  of  systems  to  be  solved  and  thus 
speeding  up  the  termination  of  the  algorithm.  We  believe  that  in  this 
way,  the  overall  efficiency  of  the  algorithm  will  be  increased. 

The  reason  we  propose  the  modified  version  of  the  algorithm  only 
for  the  case  M belonging  to  K is  in  solving  the  systems  of 

equations  = “qI  With  M € Z and  the  old 

(k) 

definition  of  the  set  Iv  , it  can  easily  be  proved  that  the  submatrix 
“n  is  Minkowski,  provided  that  the  problem  (q,M)  is  feasible;  in 
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(k+1) 


-q^  can  be  solved 


this  case,  the  system  of  equations  M^x^ 
very  efficiently  by  both  factorization  and  iterative  methods;  further- 
more, the  solution  to  the  system  of  equations  is  bound  to  be  nonnegative. 

On  the  other  hand,  with  M simply  belonging  to  Z and  our  new  definition 
(k) 

of  the  set  I'  , there  is  no  guarantee  that  the  submatrix  is  non- 
singular, as  an  easy  example  will  show,  thus  the  system  = -q 

may  have  multiple  solutions.  Although  one  can  then  solve  the  system: 

= -q  , Xjk+^  > 0,  but  it  is  not  hard  to  see  that  so  doing 
will  essentially  bring  one  back  to  the  original  Algorithm  I.  Thus 


nothing  much  is  gained.  However  if  M £ K,  then  each  system 

= -q^  is  guaranteed  to  have  a (unique)  nonnegative  solution 

(k) 

even  with  the  new  definition  of  I . Combining  the  above  facts  and 

the  advantage  (discussed  in  the  last  paragraph)  of  the  new  definition 

(k) 

of  Iv  , we  conclude  that  the  modified  algorithm  II  is  best  suited 
for  linear  complementarity  problems  with  Minkowski  matrices.  We  do 
not  recommend  it  for  problems  (q,M)  where  M € Z but  M ^ K. 

We  are  now  ready  to  state  and  justify  our  proposed  algorithm 
for  solving  the  problem  (b,c,D)  with  D £ Z and  b > 0.  The  algorithm 
is  a finite  scheme  which  requires  solving  a sequence  of  nested  sub- 
problems of  increasing  sizes.  It  starts  with  an  initialization  step 
which  determines  the  first  subproblem  to  be  solved.  An  inner  cycle 
is  then  reached  where  the  subproblem  is  actually  processed.  The  solution 
obtained  there  is  then  tested.  The  procedure  either  terminates  or 
enters  the  outer  cycle  where  a new  subproblem  (of  larger  size)  is 
determined.  The  test  for  solution  is  repeated  and  the  outer  cycle  is 
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u 


again  reached  or  else  the  procedure  terminates.  Eventually,  the 
algorithm  stops  in  a finite  number  of  steps  with  a solution  to  the 
problem. 

The  Algorithm; 

Step  0 (Initialization).  If  c^  > 0 for  every  i,  stop.  A solution 
is  given  by  x = 0.  Otherwise,  let  t = 0,  l(t)  = {i:^  < 0) 
and  J(t)  = { 1, . . . , n)\l(t) . 

Step  1 (inner  cycle).  Let  I = l(t)  and  solve  the  linear  complementarity 
problem  P(t): 


-(ci  + D^)  + DnvI  > 0 


(4.1)  v^O 


vI[-(cI  + DIlV  + DII VI]  = 0 


Let  be  "the  least  element  solution.  (See  Remark  2 below.) 


Step  2 (Outer  cycle).  If  c±  + ^(b..  " ^ > 0 for  every 

i € J(t),  stop.  A solution  is  given  by  = bj^j  - Vj.^ 

and  Xj^  = 0.  Otherwise,  let 

U(t)  = (i  £ J(t):c,  + E d,  .(b.  - v ) <0). 

1 j£l(t)  ij  0 J ~ 

Set  I(t+1)  = l(t)UIl(t),  J(t+1)  = {1,  ...,n)  \l(t+l);  replace 
t by  t+1  and  go  to  Step  1. 
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Remark  1.  This  algorithm  produces  a solution  in  solving  at  most  n 


subproblems  of  the  form  (*+.l). 


Remark  2.  Each  subproblem  P(t)  is  itself  a linear  complementarity 


problem  with  a Z-matrix,  thus  the  least  element  solution  vj^,  exists 


if  P(t)  is  feasible  (to  be  established  soon)  and  v, 


Kt) 


can  be 


obtained,  for  example,  by  Chandrasekaran' s algorithm  mentioned  earlier. 


Remark  3.  If  c < 0 in  the  original  problem,  then  the  algorithm 
reduces  to  solving  the  single  linear  complementarity  problem  (-(c  + Db),D) 
with  the  Z-matrix  D,  which  of  course  can  be  solved  very  efficiently, 
for  example,  by  Chandrasekaran' s algorithm. 


Remark  4.  The  algorithm  is  similar  to  Chandrasekaran ' s algorithm  in 
that  they  both  require  solving  subproblems  of  increasing  sizes  and  checking 
for  termination.  The  algorithms  differ  in  the  subproblems.  In  our 
algorithm,  the  subproblems  are  linear  complementarity  problems  with 
Z-matrices,  whereas  in  Chandrasekaran' s algorithm,  the  subproblems  are 
systems  of  linear  equations.  In  a forthcoming  report,  we  will  generalize 
the  problem  (b,c,D)  to  allow  some  or  all  of  the  b/s  to  be  infinity. 

This  generalization  obviously  includes  the  present  problem  and  a linear 
complementarity  pr-  jlem  with  a Z-matrix  as  particular  cases.  We  will 
also  propose  an  efficient  algorithm  for  this  generalized  problem  and 
show  that  it  unifies  the  algorithm  above  and  Chandrasekaran' s 
algorithm. 

Remark  5.  In  the  forthcoming  report,  we  will  demonstrate  that  the 
solution  generated  by  the  algorithm  is  precisely  the  least-element 
of  the  meet  semi- sublattice  S introduced  in  the  last  section. 


because  leas^  Percent  of  the  feasible  set  of  P(t+1), 

This  completes  the  inductive  step  and  the  algorithm  is  therefore 
justified. 


Monotonicity  of  the  Iterates:  The  proof  above  shows  that  the  following 

inequality  is  valid: 


'l(t+l)  - 


3H(t) 


with  the  understanding  that  the  vector  is  partitioned  in 

accordance  with  the  vector  on  the  right  side.  This  inequality  implies 
that  if  some  variable  v^  attains  the  value  0 at  some  iteration  step  t, 
then  it  remains  at  the  value  0 in  the  final  solution,  thus  can  be  dropped 
from  further  consideration.  The  subsequent  subproblems  will  then  have 
smaller  sizes,  thus  can  be  solved  more  quickly. 

In  order  to  complete  a further  analysis  of  the  algorithm,  we 
study  some  of  its  refinements  when  it  is  applied  to  solve  some  particular 
subclasses  of  problems.  We  first  consider  the  case  where  the  matrix  D 
is  tridiagonal  and  Minkowski  (i.e.  d. . = 0 if  (i-j|  > 1).  Realizations 

^ i j 

of  this  subclass  of  problems  can  be  found  in  Veinott's  taut  string 
problem  [23]  and  in  a discretized  version  of  the  unilateral  Dirichlet 
problem  with  two  obstacles  [21].  In  fact,  the  matrix  D appearing  in 
the  latter  case  has  a block  tridiagonal  structure.  This  means  that  the 
matrix  D can  be  partitioned  into  blocks  (i,  j = 1, . . . , n)  where 

D, . = 0 for  | i- j | > 1.  However,  the  diagonal  blocks  are  tridiagonal. 


Recognizing  the  storage  problems  and  difficulties  encountered  in  solving 
large  scale  linear  complementarity  problems  with  block  tridiagonal 
Minkowski  matrices  by  direct  application  of  the  algorithm,  Cottle  and 
Goheen  [6]  recently  proposed  a hybrid  algorithm  to  solve  this  class  of 
problems  which  requires  solving  subproblems  of  the  form  (*,»,Di^). 

The  algorithm  proposed  above  may  be  applied  in  this  instance.  An  investi- 
gation of  this  solution  strategy  is  reported  in  [6], 

Since  we  now  assume  that  D is  a Minkowski  matrix,  the  subproblems 
have  unique  solutions  which  must  necessarily  be  the  least  elements  of  the 
feasible  sets  of  the  respective  subproblems.  These  (unique)  solutions 
can  be  obtained  by  a number  of  efficient  algorithms.  We  r~fer  to 
Sacher[l8]  for  some  very  detailed  comparisons  of  various  methods  to 
solve  the  class  of  linear  complementarity  problems  with  tridiagonal 
Minkowski  matrices.  It  was  observed  in  the  same  reference  that  each 
principal  submatrix  of  the  matrix  D is  composed  of  s > 1 

block  diagonal  submatrices  (each  of  which  is  tridiagonal  and  Minkowski) 
with  zeroes  elsewhere.  Thus  each  subproblem  P(t)  is  again  decomposed 
into  s subproblems  each  of  which  is  a linear  complementarity  problem 
with  a tridiagonal  Minkowski  matrix  and  can  be  solved  independently  of 
the  others.  This  decomposition  of  the  subproblems  is  put  to  advantage 
in  our  adaptation  of  the  algorithm  which  is  formulated  in  flow  chart 


Notations  in  the  Flow  Chart. 


ii*  = 

11  = 

il(i)  = 1 

12  = 

13  = 


iold  = 


size  of  each  subproblem  P(t); 

index  set  identifying  the  indices  included  in  l(t); 

means  i € l(t)  and  0 otherwise; 

size  of  each  decoupled  subproblem  of  P(t); 

index  identifying  the  last  component  of  the  current 

vector  x that  has  been  determined; 

index  used  to  determine  the  current  vector  x for 

solution. 


Comments  on  the  Flow  Chart. 

It  is  composed  essentially  of  four  loops  (i)-(lV).  Loop  (I)  is 
an  initialization  of  the  indices  and  is  self-explanatory.  Loop  (II) 
consists  of  two  subloops  (ill)  and  (IV).  It  corresponds  to  an  outer 
cycle  of  the  method  where  the  subproblems  to  be  solved  are  determined. 
Each  subproblem  P(t)  is  solved  in  Loop  (ill)  and  the  current  vector 
is  checked  for  solution  in  Loop  (IV). 


Storage  Considerations. 

The  original  data  b,  c,  D are  stored  in  (5n-2)  double  precision 
numbers,  each  of  which  requires  8 bytes.  The  integer  indices  are  stored 
as  4-byte  numbers.  An  additional  n double  precision  numbers  are 
required  for  the  solution  vector  x.  Therefore  the  total  requirements, 
excluding  those  for  the  subproblems,  are  approximately 
(5n-2)  x8+nx4+nx8~  52n  bytes  of  storage.  Depending  on 


E 

[ 

(j 


J 

j 


the  method  used  to  solve  the  subproblems,  extra  storage  space  may  range 
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requirements  for  the  subproblems 


We  next  proceed  to  another  application  of  the  proposed  algorithm 


and  solve  a quadratic  program  studied  by  Cheng  [3]  in  a model  of  salary 


administration.  The  problem  is  to  find  a vector  y £ R to 


minimize  n 


subject  to  a.  < y < P 


In  the  particular  application  in  which  the  program  appeared,  the  variables 


y.  represent  employee  j’s  compensation,  while  parameters  A.  and  B 


constants  are  positive.  Making  the  substitutions,  z = A.  + B.y 


minimize  n 


subject  to  Z 


which  can  be  rewritten  in  vector  notations  as 


minimize  z Dz  subject  to  Z < z < Z 


rnmmmmmmmm 


— ..w> — 


MM'.MW'H-  , I HI-  ffii  u.,1 1 1 I I.P».|I44B,<||'1  mi  iiifVWWMFWl 


where  D=nl  -ee  € Z is  symmetric,  positive  semi -definite 
n n n 

z = (zj)>  = (Zj)  and  Z1  = (Z^) . Therefore  the  Kuhn-Tucker  conditions 

which  are  precisely  the  problem  (b, c,D)  with  b = Z1  - Z°  and 
c = TfzP , are  necessary  and  sufficient  for  global  optimality.  Note  that 
if  x solves  (b, c,D),  then  z = x + Z°  solves  (4.3). 

In  the  sequel,  we  apply  the  proposed  algorithm  to  solve  problem 
(4.3)  by  solving  the  equivalent  problem  (b, c,D).  Although  D is 
singular  (thus  D K),  each  of  its  proper  principal  submatrix  is 
Minkowski,  in  fact,  every  such  submatrix  of  order  k < n has  the  form 


“k  “ nIk  - Vk 


and  it  can  be  shown  by  an  easy  calculation  that 


1 1 /T  . 1 Tx 

Dk  n Xk  n-k  ekV 


if  k < n 


A typical  subproblem  in  the  proposed  algorithm  for  the  problem  (b,c,D) 
is  to  solve  the  linear  complementarity  problem 


- (ci + diiV  + diivi  ^ 0 


(4.4) 


Vj.  > 0 


Vj[-(ci  + Djjbj)  + D-j-jVj.]  = 0 


where  I c {1, ...,n).  Due  to  the  "simplicity"  of  the  matrix  D and 


* . * n 


30 


■ rr-~~"  "'“IWHI.M.  [HP,  Ml  ..l.i.,..  .. J, 


the  fact  that  each  of  its  proper  principal  submatrix  is  Minkowski,  we 
choose  the  modified  Chandrasekaran  algorithm  (Algorithm  II)  to  solve 
(4.4).  By  an  easy  calculation,  we  deduce  that 


c = DZ°  = nZ°  - ( E Z°.)e 
i=l  1 n 

and 


f j s Cj  + Djjbj. 

* A - 'j,  z?>vi+"<zi  - Q - Jx  (zi  - zi)em 


- nZ^  - ( E Z*  + Z Z?)e.  , , 
1 i£I  1 ljtl  1 I1' 


where  |l|  denotes  the  cardinality  of  the  set  I. 


It  is  then  necessary  to  solve  the  system  of  linear  equations 


(4.5) 


D v = f 
II. II  II  II 


where  II  is  the  subset  of  I which  consists  of  those  indices  i € I 
for  which  f^  < 0. 


*Let  lie  I,  we  denote  by  f^  the  vector  nZ^  - (Zi£I  Z^+E^j  Zi ) e | n | » 

Note  that  fT  = cT  + D b which  is  not  equal  to  c_  + D b . 

X1  il  ili  1 IlI1  Ix 

abuse  of  notation  will  occur  again  in  later  development. 


This 
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Case  1.  | II | = n.  Then  (4.5)  becomes  the  whole  system 


Dv  = f 


where  f = c + Db  . 


In  this  case,  a solution  vector  to  the  problem  (b,c,D)  is  given  by  any 


solution  to 


Dx  + c = 0,  0<x<b 


(4.6) 


D(x  + Z’”')  = 0,  0 < x < Z1  - Z°  . 


Since  the  matrix  D has  rank  n-1  and  De  = 0,  it  follows  that  such 

n 

a vector  has  a representation 


(4.7) 


x = Aefl  - Z°  where  z<^  < 4 < Z^,  i = 1,  ...,n. 


Therefore,  a solution  of  the  problem  (4.3)  is  given  z = Aen  and  the 
problem  is  thus  solved. 


Case  2.  |ll|  < n.  In  this  case,  problem  (4.5)  has  a unique  solution: 


v = D"1  f 
II  11,11  II 


*|n|  + n-jn|  e|ii|eTn 


Thus, 


(4.8) 


VI1  = ZI1 


Z zb  + Z z°. 

i£Jl  1 igl  1 

— r?im — 


e | j where  J1  = An. 


By  an  easy  calculation,  we  obtain 


(4.9) 


"fJl  + DJ1,I1VI1 


r E z + E Z 

n 71  i£Jl  igl  . 

_n  [ZJ1  n - I II I e|Jl|J  • 


We  then  look  for  non-positive  components  in  the  above  vector  and  augment 
the  index  set  II.  This  latter  step  is  repeated  until  we  arrive  at  some 
index  set  II  c I such  that  either  |ll|  = n,  in  which  case  we  have 
obtained  a solution  to  (4.3)  and  the  algorithm  stops,  or  we  obtain  a 
solution  to  problem  (4.4)  given  by  (v^*0)  where  v^  is  defined 
by  (4.8).  In  the  latter  case,  a tentative  solution  vector  to  the 
problem  (b, c,D)  is 


(4.10) 


Z zj-  + £ Z° 

„0  * 1£J1  i€J 

XI1  ' "ZI1  n - |ll|  6 | Il| 


where 


XJ1  = bJl  and  XJ  = 0 


J — (1, ..  .,n} \l  • 


Again,  by  an  easy  calculation,  we  may  deduce 

r E zS  E zl 

,0  i£J  1 i£Jl 


•V 


iWPIP  - 


The  current  vector  x is  then  checked  for  solution  by  identifying  any 
non-positive  component  in  the  vector  above.  The  algorithm  either  stops 
or  continues  with  an  augmented  index  set  I. 

Summarizing  the  above  analysis,  we  formulate  a procedure  for 
solving  problem  (4.3). 

1.  Read  in  n,  Z°  and  Z1  (N  = { 1, . . . , n} ) . 

Set  Iter (#  of  inner  loops)  = 1 and  Aver  0 = ~ Z*?. 

2.  Determine  the  set  10  = (i  £ N:Z^  < Aver  0]  and  1 10 j . 

Let  Aver  1 = £ (Z.GI0  z\  + Z^I0  Z°)  and  JO  = N\l0. 

5.  Determine  the  set  II  = {i  £ I0:Z^  > Aver  1]  and  j II | . Let 
J1  = io\ti. 

4.  If  | 11 | = n,  find  a scalar  A satisfying 

(4.11)  z°  < * < for  i = 1, . . . , n . 

In  this  case,  a solution  is  given  by  z = Aen. 

5.  If  | Il|  < n,  let  Aver  2 * (n-|ll| )_1  (Zi£J0  Z°  + £.€J1  zj). 

Determine  12  = {j  £ Jl:Z^  > Aver  2). 

J 

6.  If  12  replace  II  by  II  U 12  and  go  to  4. 

7.  If  12  = gf,  set  13  = {j  € J0:Z°  < Aver  2). 
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8.  If  13  ^ replace  10  by  10  U 13.  Update  Aver  1,  and  replace  Iter 
by  Iter  + 1.  Then  go  to  3. 


9.  If  13  = a solution  is  given  by 


2 JO  “ ZJ0’  ZJ1  ZJ1  and  ZI1  Aver  2 e|u|  • 


Remark  1.  The  procedure  shows  that  two  cases  can  occur,  namely, 

(i)  a scalar  A satisfying  condition  (4.11)  can  be  found  in  which  case, 
the  vector  having  all  components  equal  to  A is  a solution  of  problem 
(4.J),  and  (ii)  there  exist  index  sets  S (=  JO)  and  T (=  Jl)  such 
that  the  vector  z = (z^)  where  = Z°  for  i € S,  z^  = Z^  for 
i € T and  z±  = (|s|  + ItI)-1  (Zj£S  Z°  + Zj£T  zj)  for 
i € N\(S  U T)  solves  (4.3). 


Remark  2.  The  algorithm  Cheng  proposed  in  [3]  for  solving  the  quadratic 
program  (4.3)  is  essentially  a sorting  procedure  to  identify  the  sets 
S and  T if  they  exist.  A computational  comparison  of  the  two 
algorithms  will  be  reported  elsewhere. 
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5.  COMPUTATIONAL  EXPERIENCE 


This  section  is  a report  on  our  computational  experience  on 


solving  some  linear  complementarity  problems  of  the  nature  studied 


above.  Two  sets  of  experiments  were  performed;  the  first  of  which  was 
on  problems  (b,c,D)  with  tridiagonal  Minkowski  matrices11  D,  whereas 


the  second  was  on  the  quadratic  program  v.^.2)  studied  by  Cheng.  An 


was  to  test  the  capability  and 


efficiency  of  the  proposed  algorithm  in  handling  large  and  practical 


problems.  Specifically  we  wanted  to  investigate  the  number  of  sub 


problems  that  needed  to  be  solved,  how  their  sizes  grew  at  each 


iteration,  and  the  total  execution  times 


described  below,  the  modified  Chandrasekaran  algorithm  (Algorithm  II) 


was  used  for  the  subproblems 


When  complementarity  problems  of  this  kind  are  solved,  there  is 


always  the  possibility  of  fixing  some  variables  at  their  bounds  a priori 


in  order  to  reduce  the  dimensionality  of  the  problems.  This  can  be 


achieved  via  the  following  implications  which  can  easily  be  verified 


where  x is  an  optimal  vector.  We  believe  that  this  pre-processing 


procedure  will  increase  the  efficiency  of  the  algorithm 


ever,  in  the  experiments  performed  below,  this  fixing  variables  at  their 


bounds  is  not  adopted 


Applications  and  computational  results  of  the  proposed  algorithm  to  solve 
problems  having  block  tridiagonal  Minkowksi  matrices  are  reported  in  Cottle 
and  Goheen  [6], 


The  flow  chart  of  Section  4 was  coded  in  a FORTRAN  subroutine 


to  solve  problems  (b,c,D)  where 


The  computation  was  done  using  FORTRAN  H with  Opt  - 2 and  double 
precision  arithmetic  (8-bytes)  to  avoid  round-off  errors.  The  results 
of  the  first  set  of  experiments  are  summarized  in  Table  1.  The  vectors 
b and  c were  generated  according  to  the  rule 


c = -f  + gri  and  b^^  = hE^ 

where  f,  g,  h are  positive  scalars,  ri  and  °.  are  randomly 
generated  numbers  in  (0,1).  The  following  notations  are  used  in  the 
tablet 

I~q  = {i;Ci  < 0) 

I*  = {i:x*  = 0)  and  I*  = {i:xi  = b^ 

* 

where  x is  the  solution  generated.  The  sizes  of  the  subproblems 
are  listed  in  the  fourth  column.  The  size  of  the  first  subproblem  is 
always  the  same  as  |lc|  and  therefore  is  not  listed. 
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The  proposed  algorithm  has  a potential  weakness  in  handling 
problems  where  it  is  necessary  to  solve  many  subproblems  with  slowly 
increasing  size.  The  last  experiment  indicates  that  this  situation 
might  really  occur.  Nevertheless,  the  execution  (cpu)  times  in  the 
table  show  that  the  overall  performance  of  the  algorithms  is  very 
encouraging.  Further  experiments  of  the  algorithm  will  be  performed 
and  reported  elsewhere. 

The  second  set  of  experiments  is  concerned  with  the  solution  of 
Cheng's  quadratic  program  (4.2)  by  the  procedure  described  at  the  end 
of  Section  4.  Several  additional  features  of  the  procedure  sure  taken 
into  account  in  its  coding.  For  example,  the  averages  Aver  1 and  Aver  2 
can  be  updated  fairly  easily  without  computing  from  scratch.  It  is 
observed  that 


and 


Aver  1 = Aver  1 , , + _ 

new  old  n 


h ( t. 


i£I3 


Z z°) 

i£I3 


sum 


new 


sum  2 n . 
old 


E 

i£I2 


Z 


1 

i 


where  stun  2 = Aver  2 * (n  - |ll|).  The  scalar  in  step  4 caui  be 

chosen  to  be  any  number  between  max  Z*?  and  min  z}.  Finally, 

1<  i <n  l<i<n 

in  testing  12  ^ and  13  ^ $ in  steps  6 and  8 respectively,  it 

suffices  to  test  I H lnew  > I I1lold  and  1 10  I new  > I10 1 old* 

The  data  axe  generated  in  the  following  way:  for  j = 1, ...,n, 

a..  £ (0,2),  A..  € (0,2)  and  € (0,3)  are  randomly  chosen  in  the 


> 

M'  ' 


* 


IW^PPWWPpipPiPB 


intervals;  P,  = y . + Q.  where  y.  C (0,3)  is  also  randomly  produced; 
J J J * 


finally,  Z^  = A.  + B .a  and  Z^  = A.  + B,p..  The  outputs  of  the 
J J J J J «1  JJ 


experiments  are  summarized  in  Table  2.  Again,  the  cpu  times  show  very 
encouraging  results.  A comparison  of  the  performance  of  our  algorithm 
and  that  of  Cheng's  sorting  procedure  will  be  reported  elswhere. 
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D 

0 
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